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Abstract 

We suggest a method which improves the precision of studies of the primary compo- 
sition of ultra-high-energy cosmic rays. Two principal ingredients of the method are 
(1) comparison of the observed and simulated parameters for individual showers, 
without averaging over arrival directions and (2) event-by-event selection of sim- 
ulated showers by the physical observables and not by the reconstructed primary 
parameters. A detailed description of the algorithm is presented and illustrated by 
several examples. 
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1 Introduction 



The determination of the primary composition of cosmic rays with energies 
higher than ~ 10^^ eV is a real challenge. The lack of knowledge of the types of 
primary ultra-high-energy particles which induce extensive air showers makes 
it difficult to study their origin and in some cases even to determine their 
energy spectrum. More precise and less model-dependent determination of 
the cosmic-ray primary composition, especially in the highest-energy domain, 
is one of the most important tasks in contemporary astroparticle physics (for 
a review and discussions see e.g. Ref. [1]). 

Ultra-high-energy cosmic rays are observed through air showers, and to extract 
information about the primary particle, selected observables of real air show- 
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ers are compared to those of simulated ones for different primaries. Because 
of large shower-to-sliower fluctuations, one cannot determine the primary- 
particle type of an individual air shower. That is why traditional approaches 
to the composition studies are based on the determination of average char- 
acteristics of a large sample of cosmic-ray events. This approach has obvi- 
ous advantages: for a homogeneous composition, averaging smoothens out the 
fluctuations, and large statistics results in higher precision. Furthermore, the 
computational time required to perform reliable simulations of an observable 
averaged over a large sample is much smaller than that necessary for detailed 
simulations of all events in the data set. However, for mixed composition, 
averaging might become a problem because fluctuations of the discriminating 
observables over their central values are larger when air showers from different 
arrival directions are combined in one sample. Different zenith angles corre- 
spond to different atmospheric depths, and air showers detected by surface 
arrays have different stages of development and hence may have quite differ- 
ent observable parameters even for the same energies and primaries. Moreover, 
the geomagnetic field introduces the azimuthal dependence for photon [2] and 
very inclined hadron [3] showers. Even fluorescent detectors, which observe 
the shower development on its way through the atmosphere, respond to show- 
ers from different arrival directions in different ways, notably with different 
accuracies. For gamma-ray primaries at E > 5 ■ 10^^ eV, the entire shower 
development is direction-dependent. 

We suggest to improve considerably the precision of the composition studies 
by performing individual simulations for each observed high-energy event. In 
the case of a large number of events in a data sample, averaging in bins of 
arrival directions may be used. If one is not interested in the study of possible 
gamma-ray primaries with E > 5 ■ 10^^ eV, binning in zenith angle only is 
often sufficient. 

This method, in its simplest form with one observable, has been successfully 
implemented to obtain a limit on the gamma-ray fraction in the primary flux 
at E > 10^" eV using the data of the AGASA and Yakutsk experiments [4]. 
Event-by-event simulations (without selection by reconstructed parameters) 
were previously used for the composition studies in Refs. [5-7]. 

The rest of the paper is organized as follows. In Sec. 2, we present the sketch 
of the event-by-event approach to the composition studies, discuss the choice 
of the air-shower observables (Sec. 2.1), give a general idea of how to constrain 
the probable primary-particle type of an individual air shower (Sec. 2.2) and 
how to use these shower-by-shower constraints to gain information about the 
chemical composition of the primary cosmic-ray flux using a sample of events 
(Sec. 2.3). Sec. 3 contains the detailed description of the procedure outlined in 
Sec. 2 and ready-to-use formulae implementing this procedure. Sec. 4 presents 
several examples which illustrate the method by the analysis of small (and 
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hence statistically insignificant) samples of events. There, we consider only one 
composition-related observable, the muon density, and use samples of highest- 
energy AGASA and Yakutsk events as examples. In Sec. 4.1, we analyse in 
detail a single event (the highest energy air shower reported by AGASA) . The 
procedure for estimating the limit on the fraction of particular primaries in 
a given energy range is illustrated in Sec. 4.2 with the sample of six AGASA 
showers with reported energies higher than 10^° eV and known muon con- 
tent, while the algorithm to determine the best-fit composition assuming two 
possible primaries is illustrated in Sec. 4.3 with a sample of four events with 
reported energies above 1.5 ■ 10^° eV detected by the AGASA and Yakutsk 
experiments. Both samples are small and the analysis of Sec. 4 serves for il- 
lustrative purposes only. We briefly summarize and discuss novelties of our 
method in Sec. 5. Some of our notations are summarized in Appendix A. 
Appendices B and C contain technical information related to the examples 
presented in Sec. 4. 



2 Generalities 

In this section, we sketch the main elements of our approach to study individ- 
ual air-shower events and their ensembles. An operational algorithm is given 
in Sec. 3 and illustrated in Sec. 4. 

2.1 Two classes of observables 

On the basis of detector readings, several air-shower observables can be deter- 
mined experimentally; many of them are not independent of each other. For 
our purposes, we separate them into two groups which we will treat differently. 

2.1.1 E-observables 

This class contains the parameters related to the energy estimation and the 
arrival direction. For ground arrays, E-observables usually include the sig- 
nal density at a given distance from the shower core (known as 5'(600) or 
5(1000)); for fluorescent telescopes the E-observable is the total amount of 
the fluorescent light (corrected for atmospheric conditions). 

The arrival direction may be fixed because the pointing accuracy at high 
energies is often relatively high and variations of the arrival direction within 
the error bars have a negligible, compared to shower-to-shower fluctuations, 
effect on both energy estimation and measurement of the composition-related 
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parameters. In the case of poor angular resolution, for example for fluorescent 
detectors in the monocular mode, the arrival direction should be treated on 
the same footing as other E-parameters. 

In what follows, we ignore, for simplicity, the errors in the determination of 
the arrival direction and consider a single E-observable — the reconstructed 
energy of a shower, E^ec- For a fixed arrival direction, Ej-ec is in one-to-one 
correspondence with the energy estimator used by a given experiment. 

We note that E-observables are primary-dependent but this dependence is 
naturally accounted for in the method, see Sec. 2.2. 

2.1.2 C-observables 

We call the parameters used to distinguish various types of primary particles 
C-observables. Depending on the experiment, they may include muon density 
or muon richness, the slope of the lateral distribution function of the total 
signal or of muons only, shower front curvature, rise time etc. For a more 
efficient separation of different primaries, simultaneous use of several such 
parameters is justified, and we hereafter denote the set of these parameters as 

c = (ci,C2,...). 

2.2 Study of an individual event 

From the sample of air showers of interest, each event is to be studied sepa- 
rately. One simulates a number of artificial showers whose E-observables are 
consistent with the real event and which are initiated by different primaries. 
To this end, showers with different primary energies are simulated; for each of 
them, E-observables are reconstructed and compared to their observed values 
for the real shower. Simulated events are then selected for further study by 
assigning weights proportional to the estimate of how well their 13-observables 
match the data. It is important to use the same reconstruction procedure as 
in the analysis of the real data; detailed information about the detector is 
needed at this point. 

Any simplified treatment of the detector (unless quantified to be of negligible 
effect) would weaken the improvements of the proposed method over previous 
ones. 

For each simulated shower, C-observables are reconstructed, again by the same 
method as used for processing the real data. For each type of primary particle 
of interest, the distribution of C-observables of simulated showers selected by 
E-observables is obtained. 
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Finally, the C-observables measured for the real event are compared with 
these simulated distributions to determine the probabilities that the event 
was initiated by various primaries, given the measured E- and C-observables. 



2. 3 Study of an ensemble of events 

One is usually interested in a certain range of real primary energies Eq (see 
Appendix A for the summary of our notations). We will refer to this range 
as {E}\ it may be an interval {Ei < Eq < E2) or a half-line [Eq > Ei); in 
general, the restrictions on the energy may be supplemented by those on the 
arrival directions if a particular region of the sky is studied. 

One then selects the data set to be analyzed. While it should more or less cor- 
respond to the energy range of interest in terms of the reconstructed energies, 
the possibility of incorrect energy reconstruction has to be taken into account 
and the range of reconstructed energies should preferrably be extended. In an 
ideal case, all events recorded by an experiment should enter the sample, but 
most of them would contribute to the quantities of interest with almost zero 
weight. 

Each of the individual events in the real sample is to be analyzed as described 
in Sec. 2.2 but keeping track of the original energy Eq. Namely, for each event 
and each primary particle, two distributions are to be obtained for the C- 
parameters of the simulated events: (i) consistent with the observed one by 
E-parameters and having thrown energies in the domain of interest, Eq G {E}; 
and (ii) consistent with the real event by E-parameters and having thrown 
energies outside the domain of interest, Eq ^ {E}. Separate probabilities are 
to be obtained for the cases (i) and (ii). 

This is necessary because we want to constrain primary composition in the 
given domain of physical energies {E}, and contamination by showers with 
energies outside {E} should be properly taken into account. 

Then, the ensemble of these probablities obtained for all real events in the 
sample is subject to combinatoric analysis. As a result, either limits or best-fit 
parameters of the chemical composition of the primary cosmic-ray flux in the 
energy domain {E} are determined. At this last stage of the procedure one 
takes into account corrections for the "lost events". These are possible events 
with thrown energies Eq G {E} which however would escape from the sample 
either because their reconstructed energies differ strongly from Eq, or because 
of the event quality cuts. 
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3 Implementation 



In this section, we present a detailed algorithm which can be directly used in 
the analysis of the real data. 

3. 1 Study of an individual event 

Let us consider a single real event with the reported observed parameters 
(i?obs;Cobs) and the arrival direction (^,0) in horizonthal coordinates. Here- 
after, A denotes possible primary-particle type (for example, A = p refers to 
primary protons, A =Fe to iron nuclei and A = 7 to photons). For each A one 
is interested in, one generates a library of simulated showers which have: 

(i) primary particle A; 

(a) arrival direction {6, 0) (or various arrival directions consistent with {6, 0) 
within the experimental errors, as discussed in Sec. 2.1); 

(Hi) energies Eq in a relatively wide range around -Eobs- 

The thrown energies Eq of the simulated showers may be chosen randomly 
from, e.g., the interval (0.5i?obs, S-Eobs) according to I/Eq spectrum. This 
choice of the spectrum for the library enables one to control shower-to-shower 
fluctuations at high energies but saves computational time; as we will see 
below, it is not the spectrum of real particles which we assume to be real- 
ized in Nature. The interval we quote is indicative and should be adapted for 
particular events, especially at large zenith angles. 

The measurement of E-parameters is subject to statistical errors which are 
studied by experimental groups in detail. The probability distribution that 
the primary particle which produced an actual shower with the observed E- 
parameters equal to i?obs would rather produce a shower with these parameters 
equal to i^rec is denoted by gE{Ej.ac, -E'obs)- This function is usually determined 
and published by experimental groups. For instance, for the AGASA experi- 
ment QEiEj-cc, Eohs) is Gaussian in log(£'rec/-Eobs) and the standard deviation 
of E^cc is (Je ~ 0.25-Eobs [8]. Hence, in a library of simulated showers to be used 
in the analysis of a real shower with E-parameters equal to E^hs, we assign to 
each simulated shower a weight 

Wl = gE{Eohs, Ercc) ■ 

In principle, this function may depend on the type of a primary, but at this 
stage we use one and the same qe for all particles. 
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We emphasise that the function has nothing to do with physical fluctua- 
tions in the shower development, nor with primary- and direction-dependent 
systematics in the energy reconstruction, but describes instead the resolution 
power of the installation. 

Additionally, one may be interested in studying the energy spectrum different 
from the one used in the simulations of shower library (it may happen, for 
instance, that the library was simulated with the Eq"'^^^ spectrum with aub = 1 
while one is interested in say a = 2 or a = 2.7). To reproduce the required 
thrown energy spectrum, an additional weight 



is introduced (instead of E^bs in the denominator, any other typical energy 
scale may be used). 

The C-parameters are also reconstructed with some statistical errors. In ex- 
actly the same way a shower with measured C-parameters equal to c could 
produce detector readings corresponding to c'. We denote the corresponding 
probability distribution as 



Let us enumerate the showers in the library simulated for a given real event 
and for a given primary Ahj i = 1, . . . , M. The distribution of the parameters 
c for the showers consistent with the real one by E-parameters is given by 



where we denoted by wi^iA, W2,iA and CiA the values of the weights and the 
vector of C-parameters calculated for the simulated shower number i with 
the primary A, respectively. M is the normalization factor determined by the 
condition 



By definition, /^(c)(ic is the probability to have observed C-parameters in 
the region dc centered at c given the primary particle A and observed E- 
parameters. This probability assumes that a specific shower-development model 
used in the simulation code is valid. 

To compare the simulated distributions /a(c) with the observed C-parameters 
Cobs of the real event, one has to formulate a conjecture about the nature of 
its primary particle. 

One possible question is, how bad the event is described by the primary A - 
without fixing and studying other possibilities (only the shower library for A 





^c(c',c). 




(2) 
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primaries is required in this case). Given the observed values of Eohs and Cobs, 
the event is consistent with the primary A with the probabihty density 



Pa = /^(cobs)- (3) 

In the case of practical interest, when the event is unlikely being initiated by 
the primary A, the estimate of the probability that a given event could be 
initiated by the primary A is given by 

Pa^ = FA{Cobs) = J fAic)dc, k = l,2. (4) 

/A(c)</A(Cobs) 

The use of this estimate is justified if pAi -C 1. 

Another conjecture to be studied might be that the primary was either Ai or 
A2. In this case, pa^ + PA2 = 1 and 

PA, = 

Clearly, it may happen that the event is poorly described by both Ai and 
^2, that is ~ fA2 ~ and p^^ 2 are not stable. This indicates that the 
conjecture that the primary was either Ai or A2 works poorly for this event. 

In case one wishes to distinguish between several possible primaries Ai, . . . , Ak-, 
Eq. (5) should be replaced by 



Pa, 



Y.k=l /Afe(Cobs) 



3.2 Study of the ensemble of events 



Let us fix the physical energy range {E} and consider a sample of observed 
events enumerated by j = 1, . . . , iV, selected as described in Sec. 2.3. 



3.2.1 One primary type 

To constrain the fraction of the primaries A in the total flux of cosmic-ray 
particles within the energy domain {E}, one simulates libraries of A-induced 
showers described in Sec. 3.1 for each observed event in the sample. When 
determining the distribution /yi(c), one has to distinguish its two constituents, 

/A(c) = /i+^(c) + /i-)(c), (6) 
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where (c) is the distribution built with only those showers from the library 
whose thrown energies Eq belong to the domain {E}; the rest of showers 
contribute to f^\c). Consequently, one obtains for each observed event j, 
two probabilities: 

• p^^"', the probability that the event j was initiated by the primary A with 
energy Eq G {E}; in the case when the event j is unlikely to be induced by 
primary A, p^''"' is given by Eq. (4) with /a(c) = f^\c), otherwise, to get 
the conservative bound, one can set p^^'' = 1. 

• p^'*"', the probability that the event j was initiated by the primary A with 
energy £"0 ^ {E}; in the case when the event j is unlikely to be induced by 
primary A, is given by Eq. (4) with /a(c) = j\ otherwise, to get 
the conservative bound, one can set ^■^ = 0. 

Generally, these two probabilities do not sum up to unity, p^'*"' + p^^"* < 1- 
To proceed further, one needs a conjecture about other possibilities. Namely, 
we should distinguish between two other probabilities, 

• P^^'' J probability that the event j was initiated by any other primary 
than A with energy Eq G {E}; 

• P^^"'; probability that the event j was initiated by any other primary 
than A with energy Eq ^ {E}; 

Without assumptions about possible primaries other than A, the two cases 
p^^^ and cannot be distinguished by simulations. A reasonable solution is 
to assume loosely that the energy was determined correctly by an experiment, 
Eq = Ej.ec, and the probability to have the primary energy Eq follows the 
distribution gE^Eohs, Eq). Then, 

PJ^' = S- / 9E{E,^s,pE)dE, p^-)^ = s- J gE{E,^,,,E)dE, 

Ee{E} Ei{E} 

where the normalisation factor Kj is determined from the condition 

Va +Va +Pa """^a ~ 



Given the probabilities p"^^ for each real event j in the sample, one can 
determine the probability V{ni,n2) to have, among observed events, ni 
initiated by primaries A with energies Eq G {E} and n2 initiated by any other 
primaries with energies Eq G {E}. We are not interested in what happens 
outside the domain {E}, so in what follows we do not distinguish and 



p^ and use 



P — Pa "t P-^ ~ ^ Pa P-^ ■ 
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The probability V{ni,n2) can be calculated as follows. ^ The probability to 
have ii-th, . . . , 2„j-th observed events (to distinguish them, let us put them 
in order, ii < ■ ■ ■ < i„J induced by primaries A with E G {E} and ki-th, . . . , 
/c„2-th events {ki < ■ ■ ■ < A;„,2, ij ^ ki) induced by any other primaries with 
E G {E} is given simply by the product 

ij ki m„j^ij,ki 

To calculate P(ni,n2) one sums over all possible ordered subsets {{ij}, {h}) 

n<i2<--<ini 

E 

^l<^2<■■■<'^n- 

where 1 < ij,ki, rrin < A^. 

Let us suppose now that the fraction of primaries A in the total flux of particles 
with energies Eq G {E} is ca- By definition, is the probability for a single 
shower with Eq G {E} to be induced by A, while (1 — e^) is a probability for 
a single shower with Eq G {-B} to be induced by any other primary. Hence 
V{e), the probability that the observed results are reproduced for a given eA, 
can be expressed via P(ni, 77,2) by making use of the formula for a conditional 
probability, 

ni+n2<N 

neA)= E" V{n,,n,)e-^{l-eAr- (8) 

ni,n2 

(cf. Ref. [6] for a particular case rii + 722 = A^). The cases rii + 77-2 < reflect 
the possibility that N — rii — n2 events may correspond to primaries with 
Eo i {E}. 

Making use of the function P(eA), one can constrain the fraction at a given 
confidence level ^. Indeed, the allowed region of corresponds to 

P{^a) > 1 - e (9) 



^ Alternatively, a simple and fast practical way to determine V(n\,n2) makes use of 
the Monte-Carlo simulations. One generates a large number of sets of elements, 
each element marked either "A, Eq G {£;}", or "A, Eq E {£;}", or ''Eq ^ {E^ 
with the probabilities p^'*'', aiid V^~^K respectively (these probabilities are 

different for different elements j = 1,...A^). To get V{ni,n2), one simply counts 
the number of sets with rii elements marked "A, Eq G {-E}" and 77-2 elements marked 
"A, Eq G {E}" and divides this number by the total number of simulated sets. 



(7) 
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For sufficiently small values of p^'^'' , the function V{e) is a monotonically 
decreasing function and the allowed region corresponds to eA < eo? where eo 
is determined by 'P(eo) = 1 — we obtain an upper limit on eA- 

Alternatively, P{eA) may be non-monotonic, and the condition (9) may de- 
termine one or several intervals allowed for at the confidence level ^. The 
most probable value of ca corresponds to the maximal value of P{eA) over 
< ca < 1- 

Let us remind at this point that if the allowed interval for does not include 
= at a given C.L., this does not mean we can be sure that there are 
A-primaries in Nature: this latter statement is true only in the frameworks 
of the studied hypothesis (primary A versus any other primary with correct 
energy determination). 

The limit one obtains in this way is however biased because some particles 
with Eq G {E} may have reconstructed energies so different from Eq that 
the corresponding events would not enter the experimental sample chosen 
for the study. To correct the final result for these "lost" particles, another 
simulation is required. One should simulate sufficiently large number m of 
air showers initiated by the primaries A with arrival directions distributed 
randomly according to the experimental acceptance and energies randomly 
chosen from the domain Eq G {E}, this time with the realistic spectrum Eq". 
Some miost of these showers will have reconstructed energies E'rec (or some 
other parameters) such that they would escape from the sample chosen for 
the study; the fraction of these "lost" events is thus A = miost/m. The true 
fraction of A-primaries e^.tme is given by the ratio of the number of A- induced 
events M^^tme and the sum of the number of these events and the number of 
events initiated by other particles M, 

_ Ma, true 

~ Ma,,... + M ■ 

Similarly, the accessible to our study fraction of A-particles is given by 
the ratio of accessible to our study fraction of A-induced events Ma = (1 — 
A)MA,true and the sum of the number of these events and the number of events 
initiated by other particles M, 

_ Ma _ (l-A)MA,true 

Ma + M (l-A)MA,truc + M ■ ^ ^ 

Here we assume that other particles do not escape from the sample chosen for 
the study. Equations (10) and (11) enable one to place a bound on 

^^.t- = 1 - A t Ae^ ^^^^ 
from previously obtained limit on eA- 
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In the energy region to be studied, the expected spectra of different primaries 
may be different. Corresponding spectral indices may be probed within our 
approach. Indeed, in the example considered above, P(e^) depends also on 
the spectral index a^i, Eq.(l). The search for the maximum of Vie a) could 
be performed with account of this new variable thus revealing, in general, the 
excluded (or most probable, see sections below) region(s) of (e^i, a^) parameter 
space. The extension to several primaries Ai,A2,... is straightforward. 



3.2.2 Two or more primary types 

Let us turn now to the study of the conjecture that all primary particles in 
the energy domain {E} were either Ai or A2. This requires simulation and 
processing of the libraries of Ai- and A2-induced showers and no additional 
assumptions. The distributions /^'^^^^(c) are obtained for each event j and 
four probabilities are obtained: 

• p^^"* — event j initiated by the primary Ai with Eq G {E} 

• P^Ai'' — event j initiated by the primary Ai with Eq ^ {E} 

• Pa^'' — event j initiated by the primary A2 with Eq G {E} 

• p^j<^^^ — event j initiated by the primary A2 with Eq ^ {E} 



These probabilities are determined as follows. 
Now 



P^-^' = P^^' + pi? = 1 - P^Z^' - pi? (14) 
corresponds to the primaries with Eq ^ {E}. 

The probability V{ni,n2) corresponds now to the case when in the set, rii 
primaries Ai with Eq G E, n2 primaries A2 with Eq & E and N — rii — n2 
primaries Ai or A2 with Eq ^ {E} are present. It is determined by the same 
Eq. (7) with the replacements A Ai, A ^ A2, where p^ii and p*^"^-' are 
now given by Eqs. (13), (14). It also can be determined by the Monte-Carlo 
simulation described above. 

If the fraction of primaries Ai in the energy domain {E} is e^i, then, within 
the conjecture we study, the rest 1 — e^i correspond to A2. Thus Pi^Ai) is 
given by the same Eq. (8) and the region of e^^ allowed at the confidence level 
^ is again determined by P{eAi) > 1 — ^• 

The most probable value of e^i may be determined by maximization of P(eAi). 
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Account of the "lost" events requires, in this case, determination of two "lost" 
fractions, A^i and A^j. From the relation 



^ eX%i-^Aj 

e*i-(l-A^J + (l-er)(l-AAj 

one obtains 

^truc ^ — A^a) ^^^-^ 



Generalization of the procedure and equations to the case of more than two 
primary types Ak, k = 1, ... K, K > 2, is straightforward. In particular, one 
defines the probability to have a set of fractions (e^i, • • • ,^Ak)^ J^k^A^ = 1, 
among the primaries with Eq G {E}. The probability V{ni,n2) is now replaced 
by V{ni, . . . , uk) which can be calculated using a similar method. The function 
P{eAi, • • • , ^Ak) can be again expressed via V{ni, . . . , rix) by making use of 
the formulae for conditional probability, 

K K 

P{eA„...,eAK) = V{nu...,nK) -YleX , ^ e^, = 1 . 

k=l k=l 

Determination of the best-fit composition would require, in this case, maxi- 
mization of P{eAi, ■ ■ ■ ,eAii-i) with respect to A' — 1 fractions (note that 



4 Examples 



In this section, we illustrate the method with several simple (toy) examples. 
For the data, we choose the highest-energy cosmic-ray events reported by the 
AGASA [9] and Yakutsk [10] experiments (see Table B.l for the list of events 
we use and Ref. [4] for more details). The E-observable used by the both ex- 
periments is S'(600), the signal density at 600 m from the shower core. For a 
given zenith angle 6, it is in one-to-one correspondence with -Erec- Here, we 
use only one c-observable, p^(lOOO), the muon density at 1000 m from the 
core. Observed arrival directions, E^hs and p^(lOOO) are given in Table B.l 
in Appendix B, where we also present technical details of the simulation pro- 
cedure. The procedure used to reconstruct Ej-ec and p^(lOOO) is described in 
Appendix C. 

We do not discuss many experimental details and the model dependence of 
the results, so all three examples in this section should be considered as toy 
ones serving the only purpose to illustrate the method. 
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4-1 Example 1: an individual event 



The procedure described in this section has to be performed for each real 
event in the dataset. Shower hbraries for different primaries are generated as 
described in Appendix B, and ^^rec and p^(lOOO) are determined for each sim- 
ulated shower using the procedure presented in Appendix C. In this example, 
we consider the highest-energy AG AS A event (event 1 from Table B.l). 

Simulations of the AGASA detector and the reconstruction procedure [9] in- 
dicate that the function is Gaussian in \ogE with the width corresponding 
to the standard deviation of 25% of the central value in terms of energy. This 
corresponds to ai ~ 0.104 for \ogE, 

I _i£gfCE;W£obs) „+oo dE^cc 

\J iTXOi Jo £/rec 



Though the detector errors in determination of and fitting the muon LDF 
may vary with the muon density and number of detectors hit, in our toy 
example we use for the function Qc the approximation [11] of the Gaussian 
distribution with the width of 40% of the central value. In our case, there is 
only one C-observable, p^(lOOO) = c, and we thus have 

f '^2 

e 2<t2 

gc{c', c) = — -2 — , ac = 0.4c'. 

If we are interested in the primary particle of this particular event only, we 
should use Eq. (2) to determine fp and /pc- These functions are plotted in 
Fig. 1. These three functions allow to quantify, in particular, the answers to 
the questions: 

• How bad this event is described by a photon (of any energy), without fixing 
other possibilities? With the help of Eqs. (3), (4), one finds for this case 

^ ^ 0. 

• Suppose that the event may be initiated either by a photon or by a proton 
(of any energy). Which primary is prefered? Application of Eq. (5) results 
in 

p^ ^ 0.001; pp ^ 0.999. 

• Suppose that the event may be initiated either by a proton or by an iron 
nuclei (of any energy). Which primary is prefered? Application of Eq. (5) 
results in 

Pp ~ 0.69; ~ 0.31. 
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Fig. 1. Distributions of muon densities /a of simulated events consistent with 
Event 1 by the arrival direction and the reconstructed energy: thin dark line, A = 'j; 
thick gray line, A = p; dashed hne, A =Fe. The arrow indicates the observed value 
of /J^(IOOO) for Event 1. 



We emphasise once again that the probabihty Pp ~ 1 in the p-7 comparison 
does not mean it was surely a proton: this is true only for the hypothesis 
that either protons or photons can be the primaries. Indeed, Pp < 1 for the 
p-Fe comparison. The difference between protons and iron nuclei is not signif- 
icant for a single event (cf. Fig. 1). Nevertheless, this difference becomes more 
pronounced for samples of several events. 

We note that the two-primary comparison is complementary to the one-primary 
analyses because the latters assume correct energy determination for a part of 
showers (this or any other hypothesis is necessary for obtaining quantitative 
results) . 

In the following examples (Sec. 4.2, 4.3), we will study ensembles of events 
containing this Event 1. Hence we will need the functions f^'' determined in 
Eq. (6) for the energy domains {E} to be used in these examples: Eo > 10^" eV 
for Example 2 (Sec. 4.2) and Eq > 1.5 ■ 10^° eV for Example 3 (Sec. 4.3). For 
the latter case, we plot in Fig. 2 functions /^^^; the function /^"^ corresponds 
to lower energies and hence lower muon densities. The probabilities p^^^ and 
p^"* for the limit on the gamma-ray primaries, {E} = {Eq > 10^° eV}, calcu- 
lated as discussed in Sec. 3.2.1, and probabilities p^^ and p^\ {E} = {Eq > 
1.5 • 10^° eV}, calculated with the help of Eq. (13) for the proton-iron dis- 
crimination, are presented both for the Event 1 discussed here and for other 
events in the samples in Tables 1 and 3 in Sec. 4.2 and 4.3, respectively. 
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Fig. 2. Distributions of muon densities fp of simulated proton- induced events con- 
sistent with Event 1 by the arrival direction and the reconstructed energy: dotted 
Une, fp ^ (corresponds to Eq < 1.5 • 10^*^ eV); dashed hne, /p"*"^ (corresponds to 
Eo > 1.5 • 10^0 eV). 

4-2 Example 2: Limit on the gamma-ray fraction 



We concentrate here on obtaining the limit on the gamma-ray primaries with 
energies Eo > 10^° eV based on the AG AS A data. The energy domain of in- 
terest is thus {E} = {Eq > 10^° eV}. For the experimental sample we take all 
AGASA events with known muon data and energies E'obs > 8 ■ 10^^ eV: this 
widening of the energy range (cf. Sec. 2.3) compared to the domain {E} is 
justified because photon energies may be estimated incorrectly by the experi- 
ment. There are six events in the sample (events 1-6 in Table B.l); see Ref. [4] 
for the application of the same procedure to a larger dataset supplemented by 
Yakutsk events. For each of the real events we determine the probabilities p^^^ 
in the same way as described in Sec. 4.1 for the Event 1. The probabilities are 
given in Table 1. Corresponding probabilities V{ni,n2) to have ui gamma- 
ray primaries with Eq > 10^° eV and n2 other primaries with Eq > 10^^ eV, 
calculated as discussed in Sec. 3.2.1, are presented in Table 2. 

Figure 3 presents the function -P(e^) obtained from the values of V{ni,n2) 
with the help of Eq. (8). It is monotonically decreasing: this reflects the fact 
that the presence of gamma-ray primaries is disfavoured and the most probable 
fraction (corresponding to the maximal value of P{€^)) is = 0. We illustrate 
Eq. (9) for two confidence levels, ^ = 0.95 and ^ 
and < 0.47 at 95% C.L. 



0.68: < 0.21 at 68% C.L. 



To account for the "lost photons" whose true energies were Eq > 10^° eV 
but reconstructed energies were E,-ec < 8 ■ 10^^ eV, we simulated m = 1000 
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4 


0.003 


0.000 


0.887 


0.111 


5 


0.000 


0.000 


0.580 


0.420 


6 


0.000 


0.000 


0.565 


0.435 



Table 1 

Probabilities (see text for notations) for individual events 1-6 entering the 
sample for Example 2 (the limit on the gamma-ray primary fraction with {E} = 
{Eo > 1020 eV} ). 
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n2 = 1 


n2 = 2 


71-2 = 3 


n2 = 4 


n2 = 5 71-2 = 6 


ni=0 


0.000 


0.000 


0.001 


0.036 


0.234 


0.442 0.268 


ni=l 


0.000 


0.000 


0.000 


0.005 


0.009 


0.005 


ni=2 


0.000 


0.000 


0.000 


0.000 


0.000 




ni=3 


0.000 


0.000 


0.000 


0.000 






ni=4 


0.000 


0.000 


0.000 








ni=5 


0.000 


0.000 










ni=6 


0.000 













Table 2 

The probabilities V{ni,n2) to have ni gamma-ray primaries with Eq > 10^0 eV and 
712 other primaries with Eq > 10^^ eV in the sample of 6 events used in Example 2. 
Correct experimental energy determination was assumed for non-photon primaries. 
Note that ni+n2 < 6 and the rest 6 — rei — n2 particles have energies Eq < 10^^ eV 
in each case. 

photon-induced showers whose arrival directions were distributed following 
the AGASA acceptance (with the zenith angle cut of 45°) and energies £"0 > 
lO^o eV chosen randomly to follow Eq'^ spectrum suggested by several theoreti- 
cal models with raising super-GZK photon component. The fraction A ~ 0.035 
of these events had reconstructed energies E'rec < 8-10^^ eV and are therefore 
"lost". Application of Eq. (12) results in 

e^,truc < 0.22 at 68% C.L., 
e^,true < 0.50 at 95% C.L. 
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Fig. 3. P{e'y) calculated for the Example 2. Dashed line corresponds to 1 — P = 0.95 
(the 95% C.L. upper limit on e^), dotted line - to 1 - P = 0.68 (the 68% C.L. upper 
limit). 



Event 






„(+) 


PFc 


1 


0.254 


0.445 


0.136 


0.165 


2 


0.295 


0.349 


0.135 


0.221 


7 


0.163 


0.001 


0.735 


0.101 


8 


0.407 


0.107 


0.256 


0.230 



Table 3 

Probabilities Pp^^ (see text for notations) for individual events 1, 2, 7, 8 entering 
the sample for Example 3 (comparison of protons and iron nuclei). 

4-3 Example 3: Favoured hadronic composition 

To illustrate the case when two possible primary types are compared, we 
consider four events with Eohs > 1-5 ■ 10^° eV and known muon content (two 
observed by AG AS A and two by Yakutsk; events 1, 2, 7 and 8 from Table B.l). 
The primary composition in the domain {E} = {Eq > 1.5 ■ 10^'' eV} will be 
constrained within the hypothesis that the primaries are either protons or iron 
nuclei. 

In a way similar to Sec. 4.1, we process all four events and obtain probabilities 
Pp'^g, Eq. (13), which are listed in Table 3. They are transformed into prob- 
abilities V{ni,n2) to have rii proton primaries with Eq > 1.5 ■ 10^° eV and 
n2 iron primaries with Eq > 1.5 ■ 10^° eV among the four events we discuss 
(see Table 4). The function P{ep), plotted in Fig. 4, looks quite different from 
P(e^) from the previous example. There is no clear preference of either proton 
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Table 4 

The probabilities V{ni,n2) to have ni proton primaries with Eq > 1.5 • 10'^'' eV and 
n2 iron primaries with Eq > 1.5 • 10^*^ eV in the sample of 4 events used in Example 
3. Note that ni + n2 < 4 and A — ni — n2 particles have energies < 1-5 • 10^*^ eV 
in each case. 

0.5 
0.4 
0.3 

P 

0.2 
0.1 

0.2 0.4 0.6 0.8 1 

Fig. 4. P{ep) calculated for the Example 3. Notations are the same as in Fig. 3 

or iron primaries, and the function has a wide bump. The maximum corre- 
sponds to the most probable proton fraction ep ~ 0.39, but at the 95% C.L., 
the proton fraction is unconstrained. Clearly, the main reason for the weak- 
ness of the constraint is the small number of events in the sample. The limit 
at the 68% CL. illustrates the situation expected at higher confidence level 
for larger samples: the constraint reads 0.12 < < 0.73. It is reassuring that 
with only four events and for such a difficult task as to distinguish between 
different hadronic primaries, our method still allows to obtain some indicative 
information: heavier composition is slightly prefered. 

A minor additional complication in this example, as compared to the Exam- 
ple 2, relates to the calculation of the fractions of "lost" particles. The reason 
is that we unified data from two experiments in a single sample. We suppose 
that both experiments have similar sensitivities to different primaries (oth- 
erwise the unification of data sets would not be justified anyway). Then the 




95% CL. 
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expected total number of events which may be observed by an experiment 
is proportional to its exposure. We denote by ej;''"'^*-^-* and e*''"^^^^'' the proton 
fraction calculated with the help of Eq. (15) for the experiments I and II 
whose exposures are a/ and a//, respectively. The fraction ep estimated for the 
combined dataset is then 

,true(7) ^true(//) 
^true _ P |_ P 

^ l + a(//)/a/ l + ai/a(^H)' 

From the simulation of 1000 proton-induced and 1000 iron- induced showers 
from various arrival directions for each of the two cases - (I) AG AS A, 9 < 45°, 
and (II) Yakutsk, ^ < 60° - one obtains X^J^ ^ 0.05, ~ 0.06, X^J^^ ^ 

0.11, Apg^^ ~ 0.08. Since Xp and Apc are relatively close to each other, their 
contributions balance each other in Eq. (15), and the resulting limits on e*''"'' 
do not differ from the limits on ep without the correction for "lost" events 
(clearly, this would not be so in a general case). 



4-4 Testing the method with artificial samples 



To test how well the method works, we perform analyses similar to those of 
Sec. 4.1-4.3 but for a simulated data sample with known primaries. To this 
end, we generated showers with the arrival directions of events 1-8 (this al- 
lowed us to use the shower libraries generated for examples 1-3 above) and 
kept record of their primary particles and Eq. For each shower, we recon- 
structed S'(600) (and hence -Erec) and p^(lOOO) taking into account random 
detector errors. Then, we performed our analyses and compared their results 
with known primary content of the fake sample. 

The events of the simulated sample are listed in Table B.2. For the study 
similar to Sec. 4.2, we take events F1-F6 of which three are initiated by 
protons and three by photons (that is, = 0.5). Individual probabilities 
p-yjP^ for these events are given in Table 5. Application of the method gives, 
at the 95% C.L., < 0.73 if we take the sample with Eohs > 0.9 ■ 10^° eV and 
< 0.72 if we use only five events with E^hs > 10^*^ eV (the difference in the 
fractions of lost photons compensates the difference in statistics). 

For a study similar to Example 2, we take four events (Fl, F2, F7 and F8) of 
which two are initiated by protons and two by iron nuclei, that is ep = 0.5. 
The probabilities -PpFe given in Table 6. The most favoured composition 
is obtained to be ~ 0.45, but any is allowed at the 95% C.L. 
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P- 


(-) 
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Fl 


0.000 


0.000 


1.000 


0.000 


F2 


0.000 


0.000 


0.993 


0.007 


F3 


0.391 


0.287 


0.130 


0.192 


F4 


0.000 


0.000 


0.887 


0.113 


F5 


0.150 


0.534 


0.258 


0.058 


F6 


0.273 


0.090 


0.630 


0.007 


Table 5 



Probabilities p]^J^ (see text for notations) for individual events F1~F6 entering the 
artificial sample for the limit on the gamma-ray primary fraction with {E} = {Eq > 
1020 eV}. 

Event pj!^^ pf^ p'^pj p'fj 
Fl 0.281 0.084 0.450 0.184 
F2 0.159 0.421 0.086 0.333 
F7 0.159 0.064 0.300 0.476 
F8 0.458 0.070 0.323 0.149 
Table 6 

Probabilities p^^g (see text for notations) for individual events Fl, F2, F7, F8 
entering the artificial sample for comparison of protons and iron nuclei. 

5 Conclusions 

To summarise, we presented a simple method which allows to improve pre- 
cision of the studies of the primary composition of ultra-high-energy cosmic 
rays. Each event is studied individually. The simulated showers are selected 
by the physical observables to be consistent with the real event. 

Our approach exploits the fact that the uncertainty in discrimination of the 
primary-particle types by conventional methods is determined not only by in- 
trinsic fluctuations of the shower development, but also by the spread related 
to the arrival direction. Averaging over arrival directions introduces artifi- 
cial, easy-to-avoid, fluctuations in the determination of both the reconstructed 
energy (see Fig. 5 for an illustration) and composition-related parameters 
(Fig. 6). 

While the direction-related systematical effects are controlled within our method, 
there are still many other sources of uncertainty, the most important among 
them are the low precision of the measurement of composition-related quanti- 
ties in experiments and artificial fluctuations in simulated air showers due to 
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Log Eo 

Fig. 5. Direction dependence of the reconstructed energy for gamma-ray primaries. 
Plotted is the reconstructed energy (determined by the AGASA method from 
^(GOO)) versus the primary energy. Dark boxes: arrival direction of the event 1; 
crosses: arrival direction of the event 3; grey circles: arrival directions randomly 
distributed according to the AGASA acceptance (0 < ^ < 45°). Both £"0 and £Jrec 
are measured in eV. 




20.2 20.3 20.4 20.5 20.6 

Log Eo 

Fig. 6. Direction dependence of muon density for iron primaries. Plotted is p^(lOOO) 
versus the primary energy. Dark boxes: arrival direction of the event 7 (0 ~ 48°); 
crosses: arrival direction of the event 8 (0 ~ 59°); grey circles: arrival directions 
randomly distributed according to the Yakutsk acceptance {0 < 6 < 60°). Eq is 
measured in eV, is measured in m~^. 

use of the thinning approximation. Use of multiple c-observables, their care- 
ful choice and optimisation of simulations may help to improve further the 
precision of our method. 

As for any other method based on air-shower simulations, the capabilities 
of our approach are limited (and severely limited in the case of distinguish- 
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ing different hadronic primaries) because of persistent theoretical problems in 
modelling high-energy hadronic interactions. 

As compared to the traditional methods, the event-by-event study requires 
more computational time, which makes it difficult to implement it for large 
samples in the form described here. However, almost the same level of precision 
may be reached with reasonable binning in zenith (and in some cases also 
azimuth) angle. 

As we have seen in Sec. 4, the method is really strong in studies of the small 
samples, notably for the highest-energy or somewhat special (e.g. correlated 
with particular hypothetical sources) events, where statistics is insufficient to 
obtain confident conclusions with the help of traditional methods. In some 
cases, e.g. to constrain the gamma-ray primaries, even a single event may be 
succesfully studied. 
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A Notations 

A: the index enumerating different primaries, A = p, Fe, 7, . . . 
A: denotes any other primary but A. 

c = (ci,C2,...) - the set of c-observables (observables used to distinguish 
different primaries). Cobs denotes the values of these observables for a particular 
real event. 

E-observables - the observables used for the energy estimation. 
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Eq - the real energy of the primary particle; for simulated events it is the 
parameter of the simulation. 

Ej-cc - the reconstructed energy of a simulated shower; it is an E-observable 
(in the sense that it is in one-to-one correspondence with the observable used 
to estimate the energy for a given arrival direction) and can be determined 
for each air shower according to the reconstruction procedure fixed by the 
experiment (see Appendix C). In general, Ej-ec 7^ Eq due to both statistical 
and systematic errors. 

-Eobs ^ the reported energy of a given observed air shower. 

{E} - the region of primary energies Eq for which the composition is to be 
constrained. 

eA - the fraction of the primaries A in the total flux of cosmic rays with 
energies Eq G {E}. 

eA,truc ^ the fraction corrected for the "lost events" . 

A - the fraction of the "lost events", that is of the showers whose primaries 
had the energies Eq G {E} but which did not enter the sample because of 
incorrect energy determination. 



B Observed and simulated events used in examples 

Experimental information required for the simulations in Examples 1-3 of 
Sec. 4 is taken from Refs. [4,9,10] and is summarised in Table B.l. 

The artificial events F1-F8 used in examples of Sec. 4.4 have the same arrival 
directions and experimental conditions as the corresponding events 1-8 and 
are listed in Table B.2. 

For each of the eight events, we generated libraries of simulated showers in- 
duced by primaries required for the particular examples (see column (8) of 
Table B.l). For each observed event and each primary, 1000 showers were 
generated. Thrown energies Eq of the simulated showers were randomly se- 
lected between 5 x 10^^ eV and 5 x 10^° eV following the spectrum Eq^, as 
discussed in Sec. 3.1. The arrival directions of the simulated showers were 
the same as those of the corresponding real events. The simulations were per- 
formed with CORSIKA v6.204 [12], choosing QGSJET 01c [13] as high-energy 
and FLUKA 2003.1b [14] as low-energy hadronic interaction model. Electro- 
magnetic showering was implemented with EGS4 [15] incorporated into COR- 
SIKA. Possible interactions of the primary photons with the geomagnetic field 
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No. 


Exp. 


Eohs 
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</> 


p^'^'^flOOO) 


Examples 


Simulated 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


(8) 


1 


AGASA 


2.46 


36.5 


79.2 


8.9 


1,2,3 


7, V. Fe 


2 


AGASA 


2.13 


22.9 


55.5 


10.7 


2,3 


'y, v. Fe 


3 


AGASA 


1.44 


14.2 


27.5 


8.7 


2 


'y 


4 


AGASA 


1.34 


35.1 


234.9 


5.9 


2 


1 
1 


5 


AGASA 


1.05 


33.7 


291.6 


12.6 


2 


7 


6 


AGASA 


1.04 


35.6 


100.0 


9.3 


2 


7 


7 


Yakutsk 


1.60 


47.7 


180.8 


19.6 


2,3 


p, Fe 


8 


Yakutsk 


1.50 


58.7 


230.6 


11.8 


2,3 


Fe 



Table B.l 

Events analysed in examples in Sec. 4.1-4.3. The columns give (1) reference number, 
(2) experiment name (which determines the conditions for simulations), (3) the 
reported energy in units of lO^'' eV, (4) and (5) horizontal coordinates in degrees 
(azimuth angle 4> = corresponds to a particle coming from the South, (j) = 90° - 
from the West), (6) reported muon density at 1000 m from the core in units of m~^, 
(7) numbers of examples in Sec. 4 where this event is used, (8) types of primary 
particles of simulated showers. 



No. 


primary 


Eo 


-E'obs 


p°bs(iooo) 


(1) 


(2) 


(3) 


(4) 


(5) 


Fl 


P 


1.78 


2.45 


17.7 


F2 


P 


1.55 


1.84 


10.9 


F3 


7 


1.05 


0.94 


1.0 


F4 


P 


1.04 


1.35 


13.1 


F5 


7 


1.03 


1.25 


1.0 


F6 


7 


1.21 


1.77 


1.8 


F7 


Fe 


1.65 


2.10 


14.9 


F8 


Fe 


2.69 


1.71 


20.7 



Table B.2 

Artificial events analysed in examples in Sec. 4.4. The columns give (1) reference 
number, (2) primary particle, (3) the original energy in units of lO^'' eV, (4) the 
reconstructed energy in units of lO^'^ eV, and (5) reconstructed muon density at 
1000 m from the core in units of m~^. 
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were simulated with the PRESHOWER option of CORSIKA [16]. As sug- 
gested in Ref. [17], all simulations were performed with thinning level 10~^, 
maximal weight 10^ for electrons and photons, and 10^ for hadrons. 

For each simulated shower, we determined S'(600) and p^(lOOO). The plane 
orthogonal to the arrival direction has been divided into concentrical rings of 
the width of 100 m; contributions of all particles were averaged over these 
rings to reconstruct the lateral distribution function. For the calculation of 
the signal density, we used the detector response functions from Refs. [18,19]. 
The reconstructed energy - the E-observable which we compared to the exper- 
imental values - was obtained by making use of exactly the same procedure 
as used in processing the real data; see Appendix C. 



C Reconstruction of energy and muon density 

It is important to process the simulated events in a way as close as possible 
to the procedure with which the observables we use were obtained for the real 
events. For completeness, we recall here the relevant procedures for AGASA 
and Yakutsk. 

C.l AGASA energy estimation 

For the primary energy estimation, AGASA used the following procedure (the 
Takeda method, Ref. [8]). 

The lateral distribution function (LDF) of the signal is obtained by fitting the 
scincillator detectors' readings by the empirical formula [20], 

where 

r/ = 3.97- 1.79(sec^- 1), /2m = 91.6 m, i?i = 1000 m 

and r denotes the distance to the shower core. Only detectors with 300 m< r < 
1000 m are used for the fit. S'(600) is the value of the fitted LDF at r = 600 m. 

Due to the 10% one-side systematic error reported in Ref. [8], the value of 
S'(600) observed by AGASA is larger than the real one by a factor of 1.1 
which should be taken into account for the simulated events: 

5Lasa(600) = 1.15^(600) 
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The energy estimation formula is [21] 

E = 2.03 ■ 10^^ eV ^agasa(600), 

where 5*° (600) is the signal density for a vertical shower related to the real 
density of an inclined shower >S'^ga5a(600) by the attenuation formula [20] 

5Lasa(600) = 5°gasa(600) exp (sec ^ " 1) " (sec 9 - 1)^) , 

where 

Xo = 920 g/cm^ Ai = 500 g/cm^ As = 594 g/cm^ 

The simulations should be performed for the AGASA atmospheric depth, 
-^AGASA ~ 960 g/cm^ (note that Xagasa 7^ ^o)- 

C.2 AGASA muon density estimation 

AGASA detects the number of muon-counter hits, which is expected to be 
equal to the number of muons with kinetic energy i^km > 0.5 GeV/ cos 9^, 
where 6^ is the zenith angle of a particular muon. 

To determine the muon density at 1000 m from the core, the densities mea- 
sured by the detectors between 800 m and 1600 m from the core are fitted by 
the empirical LDF [22] 

where 

p = 2.52, 6 = 0.6, Ro = 266 m, Ri = 800 m. 

The value of the fitted LDF at r = 1000 m, p^(lOOO), is used as the primary 
composition estimator. 

C.3 Yakutsk energy estimation 

The energy estimation in the Yakutsk array [23,24] is quite similar to that of 
AGASA. The value of S'(600) is determined by fitting the scincillator detectors' 
readings by the following empirical formula. 
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where ^ 



i?Af = 68.0 m, i?i = 2000 m, a = 1.3, 6 = 3.24 - 2.6(1 - cos^), ^ = 3.5. 
The attenuation formula is 

5(600) = 5°(600) (^(1 -/5)exp (^-^(sec^- l)j +/3exp (^-^ (sec ^ - 1))) 
The energy conversion formula is 

E = 4.6-10i^eV p(600)J , 

where 



Xo = 1020g/cm^ A£; = 250g/cm^ = 2500 g/cm^ /3 = 0.39- [5"(600) 

The simulations are to be performed for the Yakutsk atmospheric depth, 
^Yakutsk = Xo^ 1020 g/cm^. 



C.4 Yakutsk muon density estimation 

For the two events we use, the Yakutsk muon detectors counted muons with 
Ekin > 1.0 GeV/cos^^. 

To apply one and the same procedure to both AGASA and Yakutsk results, 
we use the muon density at 1000 m as the parameter to be compared with 
simulations. The Yakutsk muon detectors have larger area and significantly 
higher saturation threshold than AGASA's, so the detector readings in a wider 
region, 400 m to 2000 m, are used in the LDF fit by the following empirical 
formula [26,23] 

/ f \ —0.75 / „ \ 0.75—6^ / „ \ —g^ 

where Rq = 280 m, Ri = 2000 m, g^^ = 8 and together with the overall 
proportionality coefficient, is a free fitting parameter. 



-0.12 



^ The actual procedure is more involved [2 5]. It takes into account the atmospheric 
conditions at the arrival time of a given event, that in particular results in an 
event dependence of Rm which varies usually within 60 m < Ri[j < 80 m. The 
detector responses then are adjusted in order to get the same energy estimate with 
Rm = 68.0 m, that allows to treat the data uniformly. Since this procedure adopted 
by Yakutsk collaboration is unambigiuos, our simple method is justified. 
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